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EFFICIENT  COMPUTATION  OF  PERIODIC  GREEN’S  FUNCTIONS 


WITH  APPLICATION  TO  GRATING  STRUCTURES 


By 

Roy  Jorgenson  and  Raj  Mittra 
University  of  Illinois 
Urbana,  IL 


ABSTRACT 


It  is  shown  that  electromagnetic  scattering  from  periodic  structures  may  be  formu¬ 
lated  in  terms  of  an  integral  equation  that  has  as  its  kernel  a  periodic  Green's  function.  The 
periodic  Green's  function  may  be  derived  from  two  points  of  view:  as  a  response  to  an 
array  of  line/ point  sources  (spatial  domain)  or  as  a  response  from  a  series  of  current  sheets 
(spectral  domain).  These  responses  are  a  Fourier  transform  pair  and  are  slowly  convergent 
summations.  The  convergence  problems  in  each  domain  arise  from  unavoidable  singulari¬ 
ties  in  the  reciprocal  domain.  A  method  is  discussed  to  overcome  the  slow  convergence  by 
using  the  Poisson  summation  formula  and  summing  in  a  combination  of  spectral  and  spatial 
domains.  A  parameter  study  is  performed  to  determine  an  optimum  way  to  weight  the 
combination  of  domains.  Simple  examples  of  scattering  from  a  one-dimensional  array  of 
strips  and  two-dimensional  array  of  plates  are  used  to  illustrate  the  concepts. 
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1.  INTRODUCTION 


The  interaction  of  electromagnetic  fields  with  periodic  structures  has  always  proven 
difficult  to  analyze.  Some  examples  of  these  structures  include  frequency  selective  surfaces 
[l],  microstrip  arrays  [2]  and  sources  inside  waveguides  [3].  One  fruitful  approach  toward 
solving  periodic  problems  involves  the  formulation  of  an  integral  equation  and  its  numeri¬ 
cal  solution  via  the  Method  of  Moments.  The  integral  equation  has  as  its  kernel  a  periodic 
Green's  function,  which  is,  unfortunately,  a  slowly  convergent  summation.  Consequently, 
the  computer  time  required  to  solve  the  problem  by  the  Method  of  Moments  is  dominated 
by  the  time  needed  to  compute  the  impedance  matrix  elements. 

In  the  past,  investigators  have  used  various  techniques  to  speed  convergence  of  the 
summation.  Functions  with  a  wide  support  in  the  spatial  domain  have  been  used  as  basis 
and  testing  functions  to  make  the  summations  in  the  spectral  domain  more  convergent 
[4.5],  Poisson's  summation  formula  [6]  has  been  used  to  speed  convergence  using  a  spatial 
domain  approach  [3.7]  or  using  a  spectral  domain  approach  [2.8.9]. 

This  report  investigates  the  efficient  computation  of  the  periodic  Green's  functions  to 
tie  together  the  methods  discussed  in  the  preceding  paragraph.  Two  examples  of  periodic 
problems  are  used  for  illustration.  The  first  example,  discussed  in  Chapter  2.  is  that  of  elec¬ 
tromagnetic  scattering  from  a  strip  grating  with  the  strips  arbitrarily  rotated  with  respect 
to  the  x  axis  as  shown  in  Figure  2.1.  The  periodic  Green's  function  for  this  example 
involves  a  single  summation.  This  simpler  problem  will  be  used  to  demonstrate  the  manner 
in  which  the  periodic  Green's  function  arises  in  the  problem  formulation,  the  slow  conver¬ 
gence  of  the  summation  and  a  method  of  accelerating  the  convergence  through  use  of 
Poisson's  summation  formula.  Results  will  be  presented  to  clarify  specific  points  of  the 
problem  and  to  show  how  optimum  parameters  are  chosen  in  order  to  maximize  computa¬ 
tional  efficiency. 


In  Chapter  3.  an  extension  of  the  techniques  developed  in  Chapter  2  is  applied  to  a 
two-dimensional  array  of  plates  on  a  skewed  coordinate  system  as  shown  in  Figure  3.1. 
The  periodic  Green's  function  arising  in  the  plate  array  problem  involves  a  double  summa¬ 
tion  and  is  computationally  more  intensive  than  the  strip  array  case.  However,  the  behavior 
of  the  two-dimensional  summation  is  similar  to  the  behavior  of  the  one-dimensional 
periodic  Green's  function.  The  conclusions  of  this  study  are  summarized  in  Chapter  4. 
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2.  THE  ONE-DIMENSIONAL  ARRAY 


2.1  Introduction 

In  this  chapter,  the  integral  equation  formulation  for  scattering  from  the  strip  grating 
shown  in  Figure  2.1  will  be  examined.  The  strips  are  rotated  to  make  an  angle  9  with 
respect  to  the  x  axis.  The  strips  are  perfect  electric  conductors  and  are  spaced  b  meters 
apart.  The  incident  field  is  a  plane  wave  where  the  direction  of  propagation  in  the  xy  plane 
makes  an  angle  of  9,  with  the  y  axis.  The  plane  wave  is  either  transverse  magnetic  (TM) 
to  z  or  transverse  electric  (TE)  to  z  . 

The  case  of  plane  wave  incidence  is  the  basis  for  analyzing  all  problems  involving 
scattering  from  a  periodic  structure.  If  an  arbitrary  field  (not  a  plane  wave)  is  incident  on 
a  periodic  structure,  no  relationship  exists  between  the  currents  of  different  unit  cells  in  the 
structure.  Therefore,  the  currents  on  the  entire  structure  must  be  treated  as  unknowns  in  a 
Moment  Method  solution.  On  the  other  hand,  if  the  incident  field  is  a  plane  wave,  then  a 
relationship  may  be  found  between  currents  of  different  unit  cells  based  on  Floquet's 
theorem:  only  the  currents  in  a  single  unit  cell  must  be  treated  as  unknowns.  This  is  dis¬ 
cussed  in  more  detail  in  Section  2.3.  The  response  due  to  an  arbitrary  source  is  found  by 
decomposing  the  arbitrary  source  into  plane  waves  and  adding  the  plane  wave  responses. 


2.2  Definitions  of  Terms 

In  the  remainder  of  this  report,  the  Fourier  transform  is  used  extensively.  The 
Fourier  transform  pair  for  the  one-dimensional  case  is  defined  as 


F(0,  )  =  /  /  (x)e  /e''dx 
/  (JC  )=  -L / F(0.  )e+jf,''dp. 


(2.1a) 


f  2. 1  b) 


where  f ( x )  is  a  function  in  the  space  domain  and  F(/ 3,  )  is  the  Fourier  transform  of  f(x) 
into  the  spectral  domain. 


<x 


The  Fourier  transform  defined  above  is  used  in  this  report  to  transform  periodic  func¬ 
tions  in  the  space  domain  into  their  equivalent  representations  in  the  spectral  domain  A 
periodic  function  may  be  viewed  as  a  convolution  of  the  function  truncated  to  one  period 
with  a  comb  function  in  space  [10],  i.e.. 

oo 

f  p(x)  =  f  (x  )*  £  S(x-mS)  (2.2) 

til  ~  — oo 

where 

/  (x  )  -  )  p  (.t  )  for  — A  <  x  <  +y 
=  0  otherwise  . 

b  is  the  period  of  Fp  (.r  )  and  *  denotes  the  convolution  operation  defined  as 

oo 

f  (x  )  *  g  (x  )  -  f  f  (x  )g  (x  —x  ‘)dx  '  (2.3) 

— CO 

Since  the  Fourier  transform  of  a  comb  function  is  also  a  comb  function,  albeit  with  a 
different  period,  and  the  transform  of  a  convolution  is  the  product  of  the  transforms,  the 
Fourier  transform  of  Equation  (2.2)  is  a  function  sampled  at  discrete  values  in  the  spectral 
domain. 

Fp(&x)  =  F(PX)  Z  8  k-^L  (2.4) 

A  concise  way  to  predict  the  location  of  the  spectral  domain  comb  components  of  Equation 
(2.4)  given  the  location  of  the  spatial  domain  comb  components  is  through  the  use  of  a 
reciprocal  lattice. 

A  periodic  geometry  has  associated  with  it  a  spatial  lattice  and  a  reciprocal  lattice  [  1 1  ]. 
The  spatial  lattice  is  a  periodic  arrangement  of  points  tn  space  and  is  formed  by  adding  an 
integer  number  (m)  of  primitive  vectors  (.S'!)  to  a  location  r  .  For  the  one-dimensional  case 
under  consideration.  S  i  =  bx  as  shown  in  Figure  2.2. 


The  overall  periodic  structure  is  formed  when  the  unit  cell,  shown  in  Figure  2.1.  is  attached 
to  each  lattice  point.  The  unit  cell  is  defined  as  the  smallest  part  of  the  structure  that, 
when  repeated,  makes  up  the  overall  structure. 


The  reciprocal  lattice  is  associated  with  the  spectral  domain  just  as  the  spatial  lattice  is 
associated  with  the  spatial  domain.  The  reciprocal  lattice  predicts  where  the  discrete  com¬ 
ponents  of  )  are  located  in  the  spectral  domain.  In  the  one-dimensional  case  under 

consideration,  the  reciprocal  lattice  is  defined  by  adding  an  integer  number  (n)  of  primitive 
reciprocal  vectors  (S  L)  to  a  location  in  the  spectral  domain  (i  ) 

K‘  =  It  +  nS,  (2.6) 

where  S ,  is  related  to  the  spatial  primitive  lattice  vector  by 

5 1  S  i  =  2rr  (2.7) 

Therefore. 


S,=  ip  (2.8) 

The  reciprocal  lattice  corresponding  to  the  spatial  lattice  of  Figure  2.2  is  shown  in  Figure 
2.3. 


The  concept  of  the  reciprocal  lattice  is  not  too  useful  for  the  one-dimensional  case. 
However,  for  two-dimensional  periodicity  on  skewed  coordinates  discussed  in  Chapter  3.  it 
provides  great  insight.  Derivation  of  the  two-dimensional  reciprocal  lattice  will  be  deferred 
until  Chapter  3. 


2.3  Formulation  of  Gp 

The  electric  field  integral  equation  will  be  used  to  solve  the  one-dimensional  strip 
array  problem.  Using  the  fact  that  the  tangential  E  field  is  zero  on  the  strip,  the  following 
equation  is  obtained: 


n  X 


—  jOJ/J.J'jis  ‘Xjp  (j  .S  ')ds‘  +  — ■ 1—  V  J  V'  J(j  Xrp  (r  .s  ’)ds  ' 


-  —  n  xE,nc  (j  )  (2.9) 


where  Gp  is  the  periodic  Green  s  function  and  J  is  the  surface  current  density  flowing  in  a 


! 
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single  unit  cell. 

The  periodic  Green's  function  arises  from  Floquet's  theorem  which  says,  given  a  plane 
wave  incident  upon  a  periodic  structure,  all  responses  will  have  the  same  periodicity  as  the 
structure  and  a  phase  shift  between  unit  cells  which  is  the  same  as  the  phase  shift  of  the 
incident  plane  wave.  For  example,  the  current  in  cell  m  is  related  to  a  corresponding  posi¬ 
tion  in  cell  0  by 

J  (r  +mbi  )  =  J  (r  )e  (2.10) 

It  is  not  necessary,  therefore,  to  consider  the  entire  structure  in  a  periodic  geometry. 

Rather,  a  single  unit  cell  may  be  considered  along  with  a  Green's  function  which  reflects  the 
relationship  of  Equation  (2.10).  The  Green  s  function  in  this  chapter  is  defined  as  the  vec¬ 
tor  potential  response  at  (x0.y0)  due  to  an  array  of  line  sources  located  at  ix'.y')  within 
each  unit  cell  and  having  a  cell-to-  cell  phase  shift  of  kx  b  due  to  the  incident  plane  wave  as 
shown  in  Figure  2.4. 


The  response  to  an  array  of  line  sources  may  be  obtained  in  two  ways.  In  the  spatial 
domain,  an  array  of  line  sources  located  at  x  '.y  '  in  each  unit  cell  may  be  represented  as 

Ja(x.y)-  Z  8(x—x'—mb)e~ji''"f’8(y—y')  (2.11) 

m  =— oo 

Summing  the  response  at  (x„.y(,)  due  to  each  line  source,  the  following  expression  is 
obtained: 


Gpixn.y^x'.y')  -  -L  £  e  ~jk '  "'h  H ,?  k  0V(x  0—x  '-mb  )2  +  (y  0~y  ’)2 
^  J  m  ——co 


(2.12) 


In  the  spectral  domain,  the  Fourier  transform  pair  is  used  to  express  the  line  array  of 
Equation  (2.11)  as  a  series  of  current  sheets.  Each  current  sheet  has  a  period  dictated  by 
the  reciprocal  lattice  and  a  cell-to-cell  phase  shift  dictated  by  the  incident  field  (Equation 
(2.6)). 


JAx.y)=  ±  Z  e  +J  e'm  ( ^8{y  —y ') 

y  m  =— oo 

where 


(2.13) 
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Adding  the  response  at  (x0.y0)  due  to  each  current  sheet,  the  spectral  domain  Green's  func¬ 
tion  is  obtained 


,  oo  o-y'1  ^  „  .  ,, 


6 

/If  3S— OO 
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(2.14a) 


where  application  of  the  radiation  condition  yields 
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(2.14b) 


To  summarize,  the  solution  of  plane-wave  scattering  from  periodic  structures  involves 
convolving  the  current  in  one  unit  cell  with  a  periodic  Green’s  function,  which  may  be 
derived  from  two  points  of  view.  Working  in  the  spatial  domain,  the  response  due  to  an 
array  of  line  sources  is  found  by  summing  the  response  due  to  each  individual  line  source 
as  shown  in  Equation  (2.12).  Alternatively,  the  spectral  domain  representation  is  obtained 
by  representing  the  array  of  line  sources  as  a  series  of  current  sheets  and  summing  the 
response  due  to  each  current  sheet  as  shown  in  Equation  (2.14).  The  two  representations  of 
the  Green's  function  are  a  Fourier  transform  pair  sampled  with  a  comb  function.  In  the  spa¬ 
tial  domain,  the  sampling  falls  on  the  spatial  lattice,  while  in  the  spectral  domain,  the  sam¬ 
pling  falls  on  the  reciprocal  lattice.  This  concept  will  be  expanded  further  in  Section  2.6. 


2.4  Convergence  Characteristics  of  Gp 

In  this  section,  the  convergence  characteristics  of  the  periodic  Green's  function  are 
examined  in  both  the  spatial  domain  (Equation  (2.12))  and  the  spectral  domain  (Equation 
(2.14)).  Using  the  asymptotic  approximation  for  the  Hankel  function  in  Equation  (2.12), 
the  spatial  domain  summation  is  found  to  behave  as 
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for  large  values  of  m.  This  expression  is  convergent  [6]  only  because  of  the  phase  variation 
in  the  numerator  of  the  summand.  For  certain  array  spacings  it  doesn't  converge  at  all  and 
for  ail  positions  of  the  basis  and  testing  functions,  it  converges  slowly. 

The  spectral  domain  formulation  of  Equation  (2.14a)  converges  rapidly  as  long  as 
y  \  In  this  report,  this  is  called  the  "off  plane"  case  since  the  observation  point  is  located 
off  the  plane  of  the  current  sheet.  The  rapid  convergence  in  the  "off  plane"  case  occurs 
because  as  m  increases,  the  plane-wave  response  to  the  current  sheets  changes  from  pro¬ 
pagating  waves  to  evanescent  waves  as  shown  in  Equation  (2.14b)  and  the  summand  decays 
exponentially.  For  the  "on  plane"  case  (y  =y ').  the  summand  no  longer  has  the  exponential 
decay  to  aid  convergence  It  behaves  as  eJ‘"x'  /m  which  converges  slowly  in  most  cases,  and 
for  some  spacings  of  the  basis  and  testing  functions  (i.e.  Ax  =0).  doesn't  converge  at  all. 

A  further  problem  in  the  spectral  formulation  is  that  for  certain  combinations  of 
array  spacing,  incidence  angle  and  summation  index  (m).  -0.  which  causes  isolated  terms 

of  the  summation  to  go  to  infinity.  Since  the  function  is  sampled  at  discrete  points,  these 
singularities  may  be  avoided  by  changing  the  angle  of  incidence  slightly.  The  individual 
terms  will  all  be  finite,  but  the  overall  behavior  of  the  function  with  respect  to  convergence 
will  not  change. 

The  reason  that  the  different  domains  exhibit  the  convergence  behavior  outlined  above 
can  be  traced  to  the  existence  of  singularities  in  each  of  the  domains.  Recall  that  the 
periodic  Green's  functions  are  a  Fourier  transform  pair  sampled  by  the  comb  function  in 
each  domain.  In  the  spectral  domain  (Equation  2.14a)  the  function  is  singular  at  the  point 
where  /3V  =  0  which  will  inevitably  occur  for  a  continuous  function  representation.  This 
implies  that  the  Fourier  transform  (the  spatial  domain  Green's  function)  is  always  a  func¬ 
tion  with  a  wide  support  and  is.  therefore,  slowly  converging.  Conversely,  the  spatial 
domain  representation  of  the  Green's  function  (Equation  (2.12))  has  a  singularity  when  the 
argument  of  the  Hankel  function  goes  to  zero.  This  singularity  is  inevitable  for  the  continu- 


ous  function  only  when  y„  =  y '  (on  plane).  For  this  case,  the  Fourier  transform  (the  spec¬ 
tral  domain  Green's  function)  is  of  wide  support  and  slowly  convergent.  As  the  (y<>— y') 
portion  of  the  argument  becomes  larger,  moving  off  plane,  the  continuous  representation  of 
the  Hankel  function  becomes  smoother  and  the  convergence  of  the  spectral  term  becomes 
more  rapid. 

To  summarize,  a  pure  spatial  domain  formulation  doesn't  converge  well  regardless  of 
the  position  of  the  basis  and  testing  functions.  The  pure  spectral  domain  formulation 
involves  a  summation  that  converges  well  as  long  as  the  basis  and  testing  functions  are  "off 
plane."  Unfortunately,  the  "on  plane"  case  inevitably  occurs.  For  example,  it  occurs  in  the 
self  term  when  the  strips  are  rotated  with  respect  to  x  (0  ^  0°  )  or  in  all  terms  when  the 
strips  are  flat  (0  =  0°). 


2.5  Smoothness  of  Basis/Test  Functions  to  Help  Convergence 

The  pure  spatial  formulation  will  be  abandoned  at  this  point  due  to  its  convergence 
problems  that  occur  regardless  of  basis/test  location.  The  pure  spectral  formulation,  which 
has  a  convergence  problem  only  in  the  "on  plane"  case,  will  be  considered  further.  It  is 
common,  in  the  pure  spectral  formulation,  to  speed  convergence  in  the  "on  plane"  case  by 
analytically  performing  the  convolution  operation  using  basis  and  testing  functions  o 
given  combined  degree  of  smoothness.  To  demonstrate  this  technique  consider  a  TM  to  i 
plane  wave  incident  on  an  array  of  flat  strips  (0  =  0°).  The  equation  for  an  element  of  the 
Method  of  Moments'  impedance  matrix  is 


,  0»»(v ') 


-dx  dx 


(2.16) 
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J  is  the  Fourier  transform  of  the  basis  function  in  the  x  direction  performed  analytically. 


T‘  is  the  complex  conjugate  of  the  Fourier  transform  of  the  test  function  also  taken  analyt¬ 
ically.  /3U„  is  defined  in  Equation  (2.13).  If  the  basis  function  is  a  pulse  located  at  the  ori¬ 
gin.  and  the  test  function  is  a  delta  function  located  at  'T  as  shown  in  Figure  2.5.  then 


)  =  Asinc  —  (2.17a) 

and 

f  '(/3vn,  )  =  e+J0'”XT  (2.17b) 

The  summand  now  behaves  as  1/m2  which  converges  regardless  of  spacing  and  quicker  than 

the  convergence  of  the  Green's  function  alone. 

Symbolically,  the  linearity  of  the  Fourier  transform  has  been  used  to  change 

T*  *  J  *  F~l  G  '2.18) 

to 


F~l  f’JG  (2.19) 

In  Equation  (2.18).  TK  (x  )  =  f(— x  )  is  needed  to  get  the  testing  function  inner  product  into 

convolutional  form.  F~l  is  the  inverse  Fourier  transform  and  takes  the  form  of  a  summa¬ 
tion  since  0,,,,  is  discrete. 

The  smoothness  of  the  basis  and  testing  functions  becomes  essential  for  convergence 
when  differential  operators  arise  in  the  integral  equation,  such  as  a  TE  to  z  plane  wave 
incident  on  an  array  of  flat  strips  (0  =  0").  In  this  case  the  left-hand  side  of  the  integral 
equation  (Equation  (2.9))  becomes 
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derivative  onto  the  test  function  converts  the  pulse  into  a  set  of  delta  functions  while  the 
transfer  of  a  derivative  onto  the  basis  function  converts  the  triangle  into  a  pulse  doublet  as 
shown  in  Figure  2.6.  Performing  the  convolution  of  the  basis,  test  and  Green’s  functions 
analytically  leads  to  the  Fourier  transform  of  a  delta  function,  which  behaves  as  1.  and  the 
Fourier  transform  of  the  pulse  doublet  which  behaves  as  1/$,  .  These  terms  together  with 
the  1//3V  behavior  of  the  Green  s  function  yield  the  same  speed  of  convergence  for  the  scalar 
potential  term  as  that  of  the  vector  potential  term  for  the  TM  case  (1/m2)  discussed  above. 
The  TE  vector  potential  term  will  converge  much  faster  (l/m4)  since  it  has  no  derivatives 
and  the  functions  to  be  convolved  are.  therefore,  smoother. 

If  the  derivatives  are  first  transferred  onto  the  Green's  function  in  Equation  (2.20), 


the  following  equation  results: 


oo  L  2  _ ft  2 

•)  T  ]dx  ‘dx 


(2.21) 


The  term  arising  from  the  vector  potential  term  behaves  similarly  to  the  TM  case  and  is 
slowly  convergent.  The  term  from  the  scalar  potential  (0,2)  doesn’t  converge  at  all.  In  this 
case,  the  smoothness  of  the  basis  and  testing  functions  is  required  to  obtain  convergence. 
Convolving  analytically.  Equation  (2.21)  becomes 


nhslrnir*'*-'1'*-' 


(2.22) 


The  level  of  smoothness  needed  for  the  above  sum  to  be  convergent  is  at  least  that  exhibited 
by  triangular  basis  and  pulse  test  functions.  This  is  the  same  level  of  smoothness  needed 
when  the  derivatives  were  transferred  to  the  basis  and  testing  functions. 

In  summary,  for  flat  strips  (0  =  0°),  the  derivatives  may  be  transferred  onto  the 
Green’s  function,  and  the  smoothness  of  the  basis  and  testing  functions  may  be  used  to  help 
convergence,  or  the  derivatives  may  be  transferred  onto  the  basis  and  testing  functions 
explicitly  and  then  the  convolution  may  be  performed.  In  either  case,  the  speed  of  conver¬ 
gence  and  the  level  of  smoothness  required  are  the  same  and  the  order  of  the  operations 
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does  not  matter. 


In  contrast  to  the  flat  case,  when  the  strips  are  rotated  (0^0"  )  with  respect  to  the  x 
axis,  the  order  of  operations  does  matter.  If  the  derivatives  are  first  transferred  to  the 
Green's  function,  the  resulting  sum  will  not  converge  regardless  of  the  level  of  smoothness 
in  the  basis  and  test  functions.  This  is  best  illustrated  by  examining  the  case  in  which  the 
strips  are  rotated  90  degrees  to  the  x  axis.  With  the  derivatives  transferred  to  the  Green  s 
function,  the  expression  for  the  the  matrix  elements  becomes 

l  >  ->  l 
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If  J  is  a  triangular  function  along  the  strip  and  T  is  a  pulse  function,  which  led  to  con¬ 
vergent  terms  in  the  flat  case,  then  for  no  overlap  between  the  basis  and  test  functions  (see 
Figure  2.7)  Equation  (2.23)  becomes 
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Note  that  the  convergence  problem  arising  from  the  scalar  potential  term  (j32)  is  mitigated 
by  the  smoothness  of  the  basis  and  testing  functions  and  the  summand  behaves  as  1/m2 
exactly  like  the  flat  case.  Also  note  that  since  this  is  an  "off  plane"  case,  the  exponential 
decav  is  the  dominant  behavior  of  the  summand. 


When  there  is  complete  overlap  between  basis  and  testing  functions  as  shown  in  Fig¬ 
ure  2.8.  Equation  (2.23)  becomes 
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The  second  and  third  groups  of  terms  in  the  braces  have  the  1/03  behavior  dictated  by  the 
level  of  basis/testing  smoothness  chosen.  Additionally,  since  these  terms  represent  the 
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contribution  of  the  basis  to  the  parts  of  the  test  not  on  the  same  plane  (labeled  1  and  2  in 
Figure  2.8).  there  is  an  exponential  decay  in  these  terms  which  also  aids  convergence.  The 
first  term  in  the  braces  does  not  converge.  It  represents  contributions  of  the  basis  functions 
to  portions  of  the  test  function  in  the  same  plane  (labeled  3  in  Figure  2.8).  Consequently, 
the  first  term  has  no  exponential  decay.  Additionally,  it  has  only  a  1/0,  behavior  rather 
than  the  1/0*  behavior  expected  from  the  level  of  smoothness  under  consideration. 

In  order  to  understand  this  problem,  the  expression 

/  .y)fJ{pxm.y  "k  ' y  'dy  dy  (2.26) 

will  be  examined  for  basis  and  testing  functions  with  various  degrees  of  smoothness,  but 

not  necessarily  the  sufficient  degree  of  smoothness  needed  to  solve  the  TE  case  under  con¬ 
sideration.  For  pulse  basis  and  delta  testing  functions  with  no  overlap.  Equation  (2.26) 
becomes 


Artnc 


-J  e*yT 


(2.27) 


Equation  (2.27)  behaves  as  1/0,  as  expected  from  a  pulse/delta  degree  of  smoothness.  The 
exponential  decay  arises  from  the  basis  and  lest  being  "off  plane"  from  one  another.  With 
complete  overlap,  the  expression  becomes 
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The  second  term  represents  the  "off  plane"  contributions  of  the  basis  to  the  test.  It  has  a 
1  / 0 v  dependence  which  stems  from  the  smoothness  of  the  basis  and  testing  functions  and 
an  exponential  decay  which  stems  from  the  "off  plane"  nature  of  the  contributions.  The 
first  term  represents  the  single  "on  plane"  contribution  from  the  point  y-0  on  the  basis. 
This  term  also  has  a  1/0,.  dependency  which  does  not  arise  from  the  smoothness  of  the 
basis  and  testing  function.  The  important  thing  to  note  is  that  the  first  term  in  Equation 
(2.28)  has  the  same  dependency  as  the  first  term  in  the  triangle/pulse  case  (Equation 
(2.25))  even  though  the  triangle/pulse  case  has  a  higher  degree  of  smoothness  than  the 
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pulse/delta  case. 

For  pulse  basis  and  pulse  testing  without  overlap.  Equation  (2.26)  becomes 
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This  expression  has  a  1/0V?  dependency  related  to  the  basis/test  smoothness  exactly  as  was 
found  in  the  flat  case.  The  exponential  decay  arises  from  the  "off  plane"  nature  of  the  prob¬ 
lem.  For  complete  overlap  the  expression  becomes 

-  -JL-sin  e~J^~  (2.30) 

j  0,  1W  2 

Again  the  second  term  has  the  expected  1/0  2  decay  from  the  smoothness  and  the  exponen¬ 
tial  decay  from  "off  plane"  contributions.  The  first  term  has  a  \/fif  decay  which  is  the  same 
as  the  "on  plane"  terms  of  the  previous  cases. 


The  pattern  that  emerges  is  as  follows:  when  there  is  no  overlap  between  basis  and 
testing  functions,  both  exponential  decay  and  degree  of  smoothness  contribute  to  rapid  con¬ 
vergence.  In  this  case  the  derivatives  may  be  transferred  to  the  Green  s  function.  For  the 
case  in  which  there  is  overlap  between  the  basis  and  test  functions  (even  touching  at  one 
point),  one  term  arises  which  behaves  as  l//3y  regardless  of  the  smoothness  of  the  basis  and 
testing  functions.  This  term  represents  the  "on  plane"  contribution  of  the  basis  to  the  test¬ 
ing  function  and.  therefore,  has  no  exponential  decay.  For  this  case,  the  derivatives  may 
not  be  transferred  to  the  Green's  function  to  obtain  a  convergent  summation.  Rather,  basis 
and  testing  functions  must  be  chosen  with  a  level  of  smoothness  to  accept  the  derivatives, 
and  the  derivatives  must  be  explicitly  transferred  onto  the  basis  and  testing  functions. 


To  ensure  that  the  above  problem  is  not  unique  to  the  0  =  90°  case,  a  strip  of  arbi¬ 
trary  rotation  9  will  be  examined  for  completely  overlapping  pulse  basis  and  pulse  test 
functions.  For  this  case,  the  "on  plane"  contribution  of  Equation  (2.26)  becomes 
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When  0  =  0"  (the  flat  case),  the  terms  in  the  braces  cancel,  but  for  all  other  cases,  the  terms 
in  the  braces  remain,  leading  to  nonconvergence  of  the  sum. 


To  summarize,  in  the  flat  case,  all  derivatives  may  be  transferred  onto  the  Green’s 
function,  and  the  smoothness  of  the  basis  and  testing  functions  ensures  the  convergence  of 
the  sum.  The  degree  of  smoothness  required  of  the  basis  and  testing  functions  is  the  same, 
whether  the  derivatives  are  transferred  onto  the  basis  and  testing  functions  or  onto  the 
Green  s  function.  A  higher  degree  of  smoothness  leads  to  faster  convergence  in  the  sum.  In 
contrast,  in  the  rotated  strip  case,  the  derivatives  must  be  transferred  onto  the  basis  and 
testing  functions  explicitly  before  convolving.  To  transfer  the  derivatives  onto  the  Green's 
f unction  leads  to  nonconvergent  terms  representing  "on  plane"  contributions  of  the  basis  to 
the  test  function.  Increasing  the  degree  of  smoothness  has  absolutely  no  effect  on  the  speed 
of  convergence.  Pulse/delta,  pulse/pulse  and  triangle/ pulse  all  have  l/f3y  terms.  Note  that 
this  is  sufficient  for  convergence  as  long  as  there  are  no  f3  terms  in  the  numerator  represent¬ 
ing  derivatives  transferred  to  the  Green's  function. 


2.6  Acceleration  of  Convergence 


A  result  of  Section  2.5  showed  that  it  is  necessary  to  transfer  all  derivatives  in  the 
scalar  potential  term  to  the  basis  and  testing  functions  when  computing  the  matrix  ele¬ 
ments  for  a  strip  grating  of  arbitrary  rotation.  This  operation  reduces  triangular  basis  and 
pulse  testing  functions  to  combinations  of  pulse  basis  and  delta  testing  functions  if  both 
basis  and  test  take  a  derivative  (Figure  2.6).  In  this  report,  the  vector  potential  term  calcu¬ 
lation  will  be  simplified  by  approximating  the  triangular  basis  by  a  pulse  with  the  same 
moment  and  approximating  the  pulse  test  by  a  delta  function  weighted  by  the  pulse  sup¬ 
port  [12]  (see  Figure  2.9).  This  approximation  can  be  justified  by  observing  that  when  the 
test  function  is  near  the  basis  function,  the  scalar  potential  term  is  the  dominant  conlribu- 
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tor  to  the  matrix  element  and  the  value  of  the  vector  potential  term  is  unimportant.  As  the 
distance  between  basis  and  testing  functions  is  increased,  making  the  vector  potential  more 
important,  the  moment  of  the  current  is  the  quantity  that  determines  the  value  of  the  vec¬ 
tor  potential.  The  moments  are  the  same  for  the  triangle  basis  and  the  approximate  pulse 
basis.  Through  the  above  approximations,  all  of  the  integrals  are  reduced  to  the  same  form: 
one  of  finding  the  response  at  the  test  location  xr  ,yT  to  an  array  of  current  pulses  of  arbi¬ 
trary  rotation  as  shown  in  Figure  2.10. 

In  order  to  accelerate  the  convergence  of  the  summation,  the  Poisson  summation  for¬ 
mula  will  be  used.  This  method  makes  use  of  the  fact  that  a  smooth,  nonsingular  function 
with  a  wide  support  in  one  domain  (either  spatial  or  spectral)  has  a  narrow  support  in  the 


other  domain.  It  also  employs  Parsevals  theorem 


fhU)f  (x)dx  =  -L fH{p)F(H)dP 


If  h(x)  is  chosen  to  be  a  comb  function  whose  elements  fall  on  the  spatial  lattice 


h  (*)=  £ 


(2.32) 


(2.33a) 


then  its  Fourier  transform  H  (/3)  is  also  a  comb  function  whose  elements  fall  on  the  recipro¬ 
cal  lattice. 


(2.33b) 


(2.34) 


H{&)=  1  £  « 

U  m  =— 03  y 

Thus,  using  Parseval  s  theorem,  a  series  may  be  represented  in  either  domain  by 

£  e~jk,mb  f  )  _  J  £  e~jL'"'t’ 8(x  —  mb  )f  (x  )dx 
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If  f(x)  has  a  wide  support  and  is  nonsingular,  implying  slow  convergence,  then  Fiji)  will 
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have  a  narrow  suppe't.  implying  rapid  convergence  and  vice  versa. 

To  demonstrate  a  fundamental  problem  with  this  approach.  Equation  (2.34)  is  applied 
to  the  slowly  convergent  "on  plane"  case  of  the  pure  spectral  domain  (Equation  (2.14a)). 


1  y  _L 

L  t  : 


2itn\ 


<  *  ’) 
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2  irm  , 
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j  e  +j  $('-*') 

V  V*o2-02  (2.35) 


In  terms  of  Parseval's  theorem. 


F(0)  = 


g  +/  >  — » ') 

2yV*.?  ~02 


(2.36a) 


1  -  +J0C 

/  (*)*  - e+^'d/3 

2wJ  2;V*,r-02 


(2.36b) 

As  discussed  previously,  the  summand  in  Equation  (2.35)  is  singular  when  012,„=*o  ■  but 
since  0,„,  is  discrete  in  m  the  singularity  is  avoidable.  In  Equation  (2.36b).  however .0  is 
continuous  and  the  singularity  cannot  be  avoided.  The  integrand  is  sharply  peaked,  so  it  is 
expected  that  application  of  Poisson  acceleration  will  not  help  convergence.  In  spite  of  this, 
if  the  integration  in  Equation  (2.36b)  is  performed. 


/  (x  )  =  -L  H  o2  ( k  u  lx0—x  '—x  I  ) 


(2.37) 


is  obtained.  Applying  Parseval's  theorem  yields 


■J—  f  f  (x  )h  (x  )dx  =  -L  £  e  jk ' mh Hn(kn\x a—x  ‘—mb  I  )  (2.38) 

This  is  the  pure  spatial  formulation  of  the  periodic  Green's  function  which  is  slowly  con¬ 
vergent.  If  the  Poisson  summation  formula  is  applied  to  Equation  (2.38).  the  result  is  the 
pure  spectral  domain  formulation  of  the  Green's  function  "on  plane."  The  unavoidable 
singularity  of  the  Hankel  function  as  the  argument  approaches  zero  leads  to  the  slow  con¬ 
vergence  of  the  "on  plane"  sum  in  the  pure  spectral  domain. 

In  both  spectral  and  spatial  domains,  application  of  the  Poisson  summation  formula 
did  not  speed  convergence  because  it  was  applied  to  a  peaked  function  with  an  unavoidable 
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singularity.  Thus,  a  better  strategy  is  to  subtract  from  the  singular  function  an  auxiliary, 
nonsingular  function  that  is  asymptotically  equal  to  the  singular  function  for  large  m;  then 
add  the  nonsingular  asymptotic  function  back  in.  The  Poisson  summation  formula  may  be 
successfully  applied  to  this  smooth,  wide,  asymptotic  function.  First,  working  with  the 
spectral  domain: 


The  first  summation  is  done  in  the  spectral  domain  and  converges  as  The  second 

summation  is  never  singular  so  the  Poisson  summation  formula  can  be  successfully  applied. 
The  second  sum  becomes 


£  e~,k'mb  Kt)(u  lx0— x  mb  I  )  (2.40) 

m  =— oo 

where  K 0  is  the  modified  Bessel  function  which  exponentially  decreases  with  increasing 
argument. 

The  operations  of  Equation  (2.18)  may  be  rewritten  using  Equation  (2.39)  as 

T*  *  J  *  F~\G)  =  TR  •  J  *  F~KG  -Ga  )  +  Ga  (2-41) 

The  inverse  Fourier  transform  F~l  is  a  summation.  The  smooth  auxiliary  function.  C“  .  has 

the  same  asymptotic  behavior  as  the  desired  function  G  and  is  summed  in  the  spatial 
domain  through  use  of  the  Poisson  summation  formula.  In  Equation  (2.41).  the  operations 
in  brackets  may  be  viewed  as  a  way  to  accelerate  a  slowly  convergent  summation  by  break¬ 
ing  it  up  into  two  rapidly  convergent  summations. 

The  convolution  operation  of  Equation  (2.41)  may  be  distributed  onto  each  domain 
and  performed  analytically  in  the  spectral  domain,  according  to 
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+  TR  *  J  *  Ga 


(2.42) 


F-'Y'J  (G-Ga)| 

This  is  equivalent  to  computing  the  impedance  matrix  elements  by  adding  the  elements  of 
two  matrices:  one  computed  in  the  spectral  domain  and  the  other  computed  in  the  spatial 
domain. 

To  accelerate  the  spatial  domain  summation  successfully  for  the  "on  plane"  case 
(y<i— y'  =  0)  the  asymptotic  behavior  of  the  Hankel  function  must  be  removed  by  moving 
the  observation  point  of  the  auxiliary  term  off  the  plane,  cb  units. 


~  £  *  A '  "'h  H  o  oV (x0~x  '-mb  )2  = 


(2.43) 


~  Y  e  A '  ”b  F  O2  \k  o'\/(x0-.r  '—mb  )2  |  -  H  „2  [fc  0V (x , x  '-mb  )2  +  cb  2 

^  J  m  =— oo 


+  Y  e  Ji ' mh  H  <?  r  oV (■*  ()— * rnA  )2  +  cb  2 

The  first  summation  remains  in  the  spatial  domain  and  converges  rapidly  because  its 
asymptotic  behavior  is  being  subtracted  out.  The  second  summation  is  brought  into  the 
spectral  domain  using  the  Poisson  summation  formula.  To  apply  Parsevals  theorem  to  the 
second  sum.  the  following  equations  are  needed: 


/  (x  )  =  -IjH ,?  J*  oV (x «)— x  '—X  )2  +  cb  2  J 

F(I3)=  r_Ltf,2Lt()Vx2  +  cb2  e-J*xdX 
J  4  j  \ 


Application  of  Parseval  s  theorem  yields 


e  '  „+,/ 3,n(>0->  ) 
2  j  (3y 


(2.44a) 

(2.44b) 

(2.45) 


where 


3  =  lULUL-k 

r*  \/n  I  *  * 


& 


* 

* 


c 


3V  = 


V*«?-02„  if*o2>0v2,„ 

-j  -JfiL-kJ  3  >  *  <? 


Since  this  is  in  the  spectral  domain  for  an  "off  plane"  observation  point,  this  sum  is  rapidly 


convergent. 


As  in  the  acceleration  of  the  spectral  domain  formulation,  the  above  procedure  may  be 
treated  as  a  way  of  quickly  summing  the  spatial  periodic  Green's  function, 

Tk*J*G=Tr*J*  (G  —  Ga  )+F~1(G'1  )  (2.46) 

By  distributing  the  convolution  operation  of  Equation  (2.46)  onto  the  different  domains  and 

performing  the  convolutions  analytically  in  the  spectral  domain,  the  following  is  obtained. 

TR  *  J  *  (G  -G°)  +  F~Kf  'J  Ga)  (2.47) 

If  cb  is  allowed  to  equal  zero,  the  asymptotic  testing  point  is  moved  "on  plane."  The  first 

term  goes  to  zero,  and  the  Green's  function  is  summed  entirely  in  the  spectral  domain. 

Since  the  test  point  of  the  asymptotic  term  is  now  "on  plane."  the  summation  is  the  slowly 

convergent  pure  spectral  domain  approach.  The  next  section  will  discuss  the  details  of 

implementing  Equations  (2.46)  and  (2.47). 

2.7  Numerical  Implementation  of  the  Spatial  Domain  Acceleration 

In  this  section,  the  details  of  implementing  the  spatial  domain  acceleration  procedure 
will  be  examined.  The  accelerated  periodic  Green's  function,  shown  in  Equation  (2.43) 
with  y o — >'  5*  0.  is  expressed  as  a  weighted  combination  of  the  spatial  domain  and  spectral 
domain. 


The  factor  that  determines  the  weighting  given  to  each  domain  is  c,  which  is  a  measure  of 
how  far  "off  plane"  the  testing  point  of  the  asymptotic  term  in  the  Green's  function  is 


located. 


Figures  2.11-2.14  show  the  nature  of  the  real  part  of  the  summand  versus  term 
number  (m)  for  various  c's.  In  this  report,  c  is  always  multiplied  by  the  cell  size  (b).  for 
example.  c=0. 1  and  cell  size  *  0.7m  moves  the  test  point  of  the  asymptotic  Green’s  function 
(cb)  0.07  m  "off  plane."  In  Figure  2.11.  c=0.001  so  the  terms  of  the  spatial  summation  are 
small  in  magnitude  and  highly  peaked  around  the  m=0  term  The  spectral  domain  has  a 
much  larger  magnitude  which  oscillates  around  zero  to  term  20  and  beyond.  For  this  value 
of  c.  the  spectral  domain  makes  most  of  the  contribution  to  the  overall  sum  and  converges 
slowly.  As  c  is  increased,  moving  further  "off  plane"  (c=0.01  and  0.1).  the  spatial  domain 
becomes  larger  in  magnitude  and  loses  its  peaked  nature.  The  spectral  domain,  on  the  other 
hand,  becomes  smaller  in  magnitude  and  more  peaked  around  m=0.  In  Figure  2.13.  for 
example,  when  c=0.1,  the  spatial  term  magnitude  is  oscillating  around  zero  until  outside  the 
m=-16:16  core  while  the  spectral  magnitude  is  zero  outside  the  m=-4:4  core.  When  c=l.  as 
shown  in  Figure  2.14.  the  spectral  domain  terms  are  essentially  zero  and  all  the  weight  is  on 
the  spatial  domain  which,  like  the  spectral  domain  of  Figure  2.11.  oscillates  around  zero  to 
term  20  and  beyond. 

The  shifting  of  weight  from  the  spectral  domain  to  the  spatial  domain  as  c  moves  "off 
plane"  is  also  seen  in  Figures  2.15-2.18.  In  this  set  of  figures,  the  value  of  the  sum  in  each 
domain  is  observed  as  the  limits  of  the  summation  are  increased  from  m*»-l:l  to  m=- 
100:100.  For  c=0.001.  the  spectral  domain  carries  all  the  weight  and  oscillates  about  its 
true  value  past  the  sum  limit  of  m=100.  As  the  asymptotic  observation  point  is  moved 
further  "off  plane"  (c=0.01.0.05  and  0.1),  the  spectral  domain  sum  converges  in  a  fewer 
number  of  terms  and  becomes  smaller  while  the  spatial  domain  sum  requires  more  terms  to 
converge  and  makes  a  larger  contribution. 

The  question  that  arises  is:  Can  the  parameter  c  be  chosen  to  minimize  the  time 
needed  to  do  the  two  summations  in  the  spatial  and  spectral  domains?  In  order  to  answer 
this  question,  a  parameter  study  was  performed  where  the  sum  limit  needed  for  conver- 


gence  to  a  given  accuracy  in  both  domains  was  plotted  versus  c  for  various  combinations  of 
cell  size,  frequency,  incident  angle  and  test  position.  What  emerged  from  this  study  is  that 
although  the  number  of  terms  needed  for  convergence  changes  with  the  parameters,  the 
general  nature  of  these  curves  remains  essentially  constant.  Two  examples  are  shown  in 
Figures  2.19  and  2.20  for  two  different  sets  of  parameters.  When  c  is  small,  the  spectral 
domain  needs  many  terms  to  converge  and  the  spatial  domain  converges  immediately.  As  c 
increases,  the  number  of  spectral  terms  needed  decreases  while  the  number  of  spatial  terms 
needed  increases  until  at  around  c=0.05.  the  graphs  cross  over  The  area  of  cross-over  is 
relatively  flat  so  c  can  be  picked  from  the  range  0.02-0.1  and  both  domains  will  be 
weighted  approximately  the  same. 

The  true  test  of  optimization,  however,  is  not  to  minimize  the  total  number  of  terms 
needed  to  perform  the  spectral  and  spatial  summations  as  was  done  above,  but  rather  to 
minimize  the  computer  time  needed  to  perform  the  calculation  in  Equation  (2.46)  or  Equa¬ 
tion  (2.47)  applied  to  the  geometry  of  Figure  2.10.  Prior  to  examining  these  results.  Equa¬ 
tions  (2.46)  and  (2.47)  must  be  discussed  in  greater  detail. 

The  implementation  of  Equation  (2.46)  is  subsequently  called  Method  1.  A  numerical 
Romberg  integration  routine  is  used  to  integrate  the  prime  coordinates  of  Gp  over  the  one¬ 
dimensional  pulse  in  the  unit  cell.  Since  the  test  is  a  delta  function,  the  convolution  with 
the  test  function  becomes  an  evaluation  of  J*  [G^,  ]  at  the  point  xr  ,yr .  For  each  s  '  chosen  by 
the  integration  routine  along  the  strip,  the  spectral  and  spatial  summations  of  Equation 
(2.46)  are  summed  to  accuracy.  First  the  core  (m— 2:2)  is  summed  in  each  domain  to  deter¬ 
mine  which  domain  is  the  dominant  contributor  to  the  Green's  function.  The  dominant 
contributor  establishes  an  absolute  accuracy  of  the  summations  to  minimize  the  time  spent 
in  computing  a  summation  that  has  an  insignificant  contribution  to  the  integrand.  When 
the  test  is  coincident  with  the  basis  function,  the  singularity  is  removed  from  the 
H (k  „  \  x  —x  '—mb  I)  term  and  computed  analytically.  The  singularity  does  not  occur  in 
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the  asymptotic  terms  since  c>0. 

The  implementation  of  Equation  (2.47)  is  subsequently  termed  Method  2.  In  this 
method,  the  basis  and  testing  functions  are  distributed  onto  the  spatial  and  spectral  domain. 


fTU)fjU-)±  1 

4  J  />>  88- 


(2.49) 


H  o  j*  „V(x  0—x  '—mb  )2  +  (y  0— y  *)2  —  ■//«  jfcoVCxo- ■*  mb  )2  +  ( ly„— y  '  I  +cb  )2  j  ds  'ds 


-j  < 1  y  o->' ' 1  +cft 


-dy  'dy 


The  spatial  domain  integrals  are  done  numerically,  as  in  Method  1.  while  the  spectral 
domain  integrals  are  performed  analytically.  When  no  overlap  exists  between  J  and  T  in 
y  (see  Figure  2.21).  the  spectral  domain  sum  becomes 


J  oc  J  sin(0„„  cosO ±  fiy  sin0)— 
b  ,n  J  P?  (/3U„  cosfl±0,  sin0) 


,+/  'T~&\  l  vr  '  )  „  < 


yr  >0 


(2.50) 


When  J  and  T  do  overlap  in  y.  the  spectral  domain  summation  becomes 


1  V  _+,*v>r  .+""‘'(2.51) 

jA  i  jA  2  2  jfiy 

A  i  =  0  cosO— 0y  sinO 
A  2  -  0,,„  cos 0  +  Py  sin0 

Note  that  when  overlap  exists  between  the  basis  and  test  functions,  certain  terms  of  Equa¬ 
tion  (2.51)  decay  as  1/0, r  and  have  no  exponential  decay.  In  order  to  obtain  exponential 
convergence,  therefore,  it  is  necessary  to  move  "off  plane"  enough  so  that  no  overlap  occurs 
between  the  basis  and  test  functions.  In  Method  1.  the  asymptotic  test  point 
(  I  v , , — y  ’  I  +cb  )  was  redefined  for  every  point  called  by  the  integration  routine  because  the 
spectral  and  spatial  contributions  were  integrated  together.  In  Method  2.  the  asymptotic 
test  point  must  be  fixed  for  the  entire  calculation  because  the  spectral  and  spatial  contribu- 


tions  are  integrated  separately.  This  means  that  in  Method  2.  it  is  possible  for  the  asymp¬ 
totic  terms  to  be  singular  if  c  is  chosen  such  that  V (x , x  '—mb  )2+(  I  v , y  '  I  +cb  ) 2  falls  on 
the  basis  function  as  shown  in  Figure  2.21. 


Method  3  is  similar  to  Method  2  except  that  the  current  pulse  in  the  spatial  domain 
integral  outside  a  core  region  (m*-l:l)  is  approximated  as  a  delta  function  weighted  by  the 
support  of  the  pulse.  Inside  the  core,  where  the  integrand  varies  quickly,  the  integration  is 
performed  numerically.  Equation  (2.49)  becomes 


fns,yfjw)±  Z  (2.52) 

*  J  m  =-l 

H |&,>>/ (x  —x  '—mb  )2  +  (v  — y  ‘)2  J  —  H  2  |&(>%/ (x  —x  ’—mb  )2  +  (  ly  — y  ’  I  +cb  )2 

+  A  £  e~jl'mh 
^  j  m  *-1.0+1 


H 0  oVC.t*  —x,  —mb  )2  +  (y,  — y,  )2  —  Hi  k  ()V (xt  —x,  —mb  )2  +  (  I  yk  —y,  I  +cb  )2 


^  m  oo  **  J  P  y 


Figures  2.22-2.26  show  the  time  needed  to  compute  the  convolution  of  basis,  test  and 
periodic  Green's  functions  using  Methods  1.2  and  3  when  the  strips  are  flat  (0  =  0°).  As  in 
Figures  2.19  and  2.20.  although  the  calculation  time  changes  with  parameters  such  as  test 
location,  frequency,  array  spacing  and  incident  angle,  the  shape  of  the  curves  remains  essen¬ 
tially  the  same.  Method  3  is  the  fastest  method  regardless  of  the  parameters,  but  since  it 
involves  an  approximation  in  the  spatial  domain,  it  is  not  as  accurate  as  Methods  1  and  2. 
This  inaccuracy  becomes  more  pronounced  as  c  increases  and  the  spatial  domain  gets  more 
weight.  Method  1  shows  a  shape  predicted  by  Figures  2.19  and  2.20.  If  c  is  too  small 
(c<0.01).  too  much  time  is  spent  summing  the  spectral  domain  and  the  required  time  for 
the  calculation  increases.  As  c  increases  (0.01  <c<0.08).  the  time  goes  to  a  minimum  then 
slowly  increases  as  the  spatial  domain  becomes  over-weighted. 
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Figure  2.22  Time  required  for  Methods  1.2  and  3  vs.  c  (flat  case) 
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Figure  2.23  Time  required  for  Methods  1.2  and  3  vs.  c  (flat  case) 
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Figure  2.24  Time  required  for  Methods  1.2  and  3  vs.  c  (flat  case) 


15.  OOC- 


The  time  requirements  for  Method  2  are  similar  to  the  requirements  for  Method  3  for 
small  c  but  increase  with  increasing  c.  crossing  over  Method  1  at  around  c-0.02.  Method  2 
time  is  dominated  by  the  time  needed  to  integrate  numerically  in  the  spatial  domain. 
Although  Method  2  has  a  simpler  integrand  than  Method  1,  it  is  difficult  to  specify  an  abso¬ 
lute  accuracy  correctly  for  the  summation  in  Method  2.  The  accuracy  is  overspecified  in 
Method  2  and.  therefore,  requires  more  time  than  Method  1.  In  Figure  2.26.  at  c-0.077. 
Method  2  exhibits  a  time  spike.  This  occurs  because  the  real  part  of  the  spatial  integral  is 
approximately  zero,  and  the  Romberg  integration  routine  calls  the  integrand  many  times  in 
order  to  get  a  sufficient  degree  of  relative  accuracy.  As  seen  in  Figures  2.27  and  2.28.  the 
integrand  itself  is  very  smooth  and  the  large  number  of  calls  is  unnecessary.  An  absolute 
accuracy  parameter  could  be  specified  in  the  integration  routine  to  alleviate  this  problem, 
but  this  was  not  done. 

The  self-term,  shown  in  Figure  2.29.  deserves  special  consideration.  Methods  2  and  3 
exhibit  the  same  behavior  as  before.  Method  1  takes  far  more  time  for  c<0.01  than  could 
be  explained  by  saying  that  the  spectral  sum  is  overweighted.  The  explanation  for  this 
behavior  comes  from  an  examination  of  the  integrand.  When  c  is  close  to  the  strip,  the 
integrand  is  ill-behaved,  as  shown  in  Figure  2.30.  The  singularity  has  been  subtracted  only 
from  the  nonasymptotic  term  in  the  Green’s  function.  When  c  is  small,  however,  the 
asymptotic  terms  are  also  tending  to  be  singular.  Moving  "off  plane"  a  bit  more,  as  shown 
in  Figure  2.31.  causes  the  integrand  to  become  better  behaved. 

When  the  strip  is  rotated  (9  =  45°  ).  the  sum  has  essentially  the  same  behavior  as  in 
the  flat  case  with  the  exception  of  two  features  (see  Figure  2.32).  The  first  feature  is  that 
Method  1  no  longer  increases  in  time  when  c<0.01  because  since  the  strip  is  rotated,  most 
of  the  points  called  by  the  integration  routine  are  farther  "off  plane"  than  the  specified  "off 
plane"  factor.  The  second  feature  is  the  drop  in  time  exhibited  by  both  Methods  2  and  3  at 
c=0.()5.  This  occurs  because  c  has  moved  from  "on  plane."  where  the  spectral  convergence 


Figure  2.27  Real  pari  of  integrand  for  time  spike  in  Method  2.  Figure  2.26 
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Figure  2.28  Imaginary  part  of  integrand  for  time  spike  in  Method  2.  Figure  2.26 
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behaves  as  1/0'  to  "off  plane"  (c  x.7>  OSsin  (45°  ))  where  the  convergence  is  exponential. 

As  c  approaches  zero  in  Method  2.  the  accuracy  of  the  sum  must  be  specified  more  pre¬ 
cisely  because  the  sum  is  behaving  as  1/m2  with  no  exponential  decay.  To  study  the  region 
when  c  approaches  zero.  Method  2  was  summed  until  four  digits  of  accuracy  were  obtained 
in  the  final  answer  for  all  c.  The  result  is  shown  for  a  flat  strip  (0  =  0°)  case  in  Figure  2.33 
and  for  a  rotated  (0  =  45° )  case  in  Figure  2.34.  For  the  flat  case,  the  best  choice  for  c  was 
found  to  be  c=0  for  the  given  size  of  basis  and  testing  functions.  No  weighting  in  the  spa¬ 
tial  domain  is  necessary  for  Method  2.  because  the  smoothness  of  the  basis  and  testing  func¬ 
tions  help  convergence  for  all  combinations  of  these  functions.  Since  the  convolution  is 
done  analytically  in  the  spectral  domain,  there  is  no  numerical  integration  involved.  When 
c  ^0.  a  numerical  integration  must  be  performed  which  dominates  the  calculation  in  time 
even  though  the  contribution  from  the  integration  is  small.  In  the  rotated  case,  the  best 
choice  for  c  is  0.05  <c  <0.15.  Here,  smoothness  of  basis  and  testing  functions  does  not  help 
convergence  in  the  spectral  domain.  In  order  to  get  exponential  convergence,  we  must  go  "off 
plane"  (c>0.05).  In  this  case  the  time  needed  for  numerical  integration  does  not  outweigh 
the  time  needed  to  sum  in  the  spectral  domain  accurately. 
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Figure  2.34  Time  required  for  Method  2  vs.  c  for  4  digit  accuracy  (rotated  case) 


3.  THE  TWO-DIMENSIONAL  ARRAY 


3  1  Introduction 

In  this  chapter,  the  concepts  developed  in  Chapter  2  will  be  extended  to  examine  the 
formulation  of  scattering  from  the  two-dimensional  array  of  plates  shown  in  Figure  3.1. 
The  plates  are  aligned  perpendicularly  to  the  xy  plane  and  rotated  to  make  an  angle  9  with 
respect  to  the  x  axis.  The  plates  are  arranged  along  a  skewed  axis  S  j  and  S2 ■  The  incident 
field  is  a  plane  wave  with  a  direction  of  propagation  9,  with  respect  to  z  and  <f>,  with 
respect  to  x 

3  2  Definition  of  Terms 

The  Fourier  transform  needed  for  the  two-dimension  array  is 

f'U 3(  ,0y  )  =  JJf  (x  ,v  '  +<5'r  'dxdy  (3.1a) 

/  •*  >  -  -rA t/ f F  <0,  .0,  )«*'<».  •  *».  >  <d  d  (3.1b) 

v  Itt  r 

where  f(x.v)  is  the  function  in  the  space  domain  and  F(0(  .(3V  )  is  the  Fourier  transform  of 
f(x.v)  into  the  spectral  domain. 

The  spatial  lattice  for  this  problem  is  shown  in  Figure  3.2  is  defined  by  use  of  a  trans¬ 
lation  vector  p,„„ 

=  P  +  P,„n  (3.2) 

=  p  +  mS !  +  nS  2 

w  here  S  1  and  S2  are  the  primitive  vectors  defined  as 

S !  =  cv  (3.3a) 

S2  —  dcosClx  +  d  sin  fly  (3.3b) 

Therefore,  the  translation  vector  in  Cartesian  coordinates  is 

P„m  =  (nd  cos (l  )x  +  (me  +  nd  sin fl  )v  (3.4) 

The  reciprocal  lattice  is  defined  through  the  use  of  a  reciprocal  translation  vector  Kmn  . 
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K=  K  +  Kmn 

—  K  +  mS  i  +  nS  2 

The  reciprocal  primitive  lattice  vectors.  5  j  and  S2.  are  defined  such  that 

Sj  5,  =  2 n  S;S,  =  0 
S  i  S2  =  0  S2  S2  =  2ir 


(3.5) 


(3  6) 


Therefore. 


S  i  =  — ^Tr~--  ly  cosfl  —  x  sinfl 
ccostl  r 


2ir 


d  cosfl 

and  the  reciprocal  translation  vector  in  Cartesian  coordinates  becomes 


(3  7a) 

(3  7b) 


=  2v 


n  _  m  sinfl 
d  cosf)  ccosfi 


•  .  2  lrm  - 

X  +  - V 


(3  8) 


The  reciprocal  primitive  lattice  vectors  are  shown  in  relation  to  the  spatial  primitive  lattice 
vectors  in  f  igure  3  3 

3  3  Formulation  of  (Jr 

The  electric  field  integral  equation  ( Fquaiion  (2  6 ) )  will  be  used  to  analvze  the  doubly 
periodic  array  of  plates  In  this  case  the  Green  s  (unction  is  defined  as  the  vector  potential 
response  to  an  array  of  point  sources  In  the  spatial  domain  an  array  ol  point  sources 
located  at  x  v  '  .z  in  each  unit  cell  may  be  represented  as 


JJ  t  .y  JT  )  =  £  £  6(W-P„„  k  ■*  ‘"Hz  ) 

n  —  — >*>n  — 

T  he  response  at  r,,.v,,.c,,  to  each  point  source  mac  be  summed  to  obtain 


i  3  «>) 


C,  <  r„.r-)*  £  £ 


•i  'J  1  n  ^ 


(  1  MM 


i  *  — -»«et  —  - 


In  the  spectral  domain,  a  point  source  arrav  may  be  expressed  as  a  double  summation 
)!  surrent  sheets  through  the  use  ol  the  Fourier  iranslorm  pair  '1  quation  '  3  I  >)  Fach  >1 
the  current  sheets  has  a  period  dictated  bv  the  reciprocal  lattice  and  a  cell  i<  vel I  phase 


v 
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shift  dictated  by  the  incident  field. 


Ja  U  .y.z)=-L-  t  t  e*J  (*mn'r‘nc  )<?_?)S(z  -z  ') 


C  A 

v  ■  tn  =— oo  n  =— oo 


(3.11) 


C.A.  is  the  area  of  the  unit  cell.  Kmn  is  defined  in  Equation  (3.8)  above  and  kmc  is  the  pro¬ 
pagation  constant  of  the  incident  wave.  Adding  the  response  at  x0,y„,z0  of  each  current 
sheet  the  following  equation  is  obtained: 


i  °°  00  _  z  tIIK)(p0  p) 

c'(r',ir)=cX  2  Z  - - - 

/))  =—00 n  s=— 00  J  7 


(3.12) 


where  application  of  the  radiation  condition  yields 


V*  ,r  -dr-  0:-  if  k  2  >  0  2  +  0? 
y~  -jyftfTJFkJ  if  0,2  +  0/>*,? 


(3.13) 


a  /  rj  p  \  *  •%  fH  Sin  O  1 

0,  =  (K„m  -k,nc  )  x  -  2tt  dcQs(i  -  ccQsn 

0y  =(Kmn-k,nc)y  =  — 

c 

The  spatial  domain  formulation  of  Gp  converges  slowly  as  explained  in  Chapter  2. 
The  spectral  formulation  converges  rapidly  when  z„^z'  (the  "off  plane"  case)  and  con¬ 
verges  slowly  when  z,,=z  ’  (the  "on  plane"  case).  As  with  the  strip  array,  since  the  plates  are 
not  flat  on  the  xv  plane,  the  derivatives  of  the  scalar  potential  term  may  not  be  transferred 
onto  the  Green  s  function,  rather  the  derivatives  must  be  transferred  explicitly  onto  the 
basis  and  testing  functions.  In  order  for  the  subsequent  integrations  to  make  sense,  rooftop 
basis  and  razor  testing  function  order  of  smoothness  is  chosen  as  shown  in  Figure  3.4.  This 
is  the  three-dimensional  analogue  to  the  triangle  basis  and  pulse  testing  functions  used  in 
Chapter  2.  If  both  basis  and  test  functions  take  a  derivative,  the  rooftop/ razor  reduces  to 
combinations  of  two-dimensional  pulses  and  delta  testing.  For  the  vector  potential,  the 
rooftop  is  approximated  as  a  two-dimensional  pulse  with  the  same  moment  as  the  rooftop 
and  the  razor  as  a  delta  function  weighted  by  the  support  of  the  razor.  All  integrals  are. 
therefore,  reduced  to  the  same  form:  the  response  at  a  point  xT  ,yT  .:T  due  to  an  array  of 


two-dimensional  current  pulses  located  on  the  spatial  lattice. 


3.4  Acceleration  of  Convergence  in  Spatial  Domain 

In  order  to  accelerate  the  convergence  of  the  spatial  domain  sum  (Equation  (3.10)),  the 
asymptotic  behavior  of  e~jip  /R  is  added  to  and  subtracted  from  the  periodic  Green  s  func¬ 
tion  by  moving  off  the  xy  plane  cCA  units. 


co  oo  0  V  *  *  +  (z  )  oo  qo  r* 

I  Z  P-"_ U - -  =  Z  I 

4t ry  I  Pu—p'—p,,,,,  I  2  +  (.2  0—z  'Y  m  =-<*> 


3 — oo/i  s-ai 


(3.14) 
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it*  n  ^ 


-Jl  o  v  1  o0-p~pmn  »  2  ♦  '  1 2  0-2  1 )  +cC A  )2 


+  E  I 


n  =-oo/i  =— x>  47r"v  1  Pi>~P  ~P„„  1 3  +  (  I  c :  ' !  +cC/l  )2 

The  first  summation  remains  in  the  spatial  domain  and  converges  rapidly  because  the 


asymptotic  behavior  is  subtracted  out.  The  second  summation  is  smooth,  nonsingular  and 
slowly  converging  It  is  brought  into  the  spectral  domain  by  means  of  Poisson  summation 
formula 

l  or  a  two-dimensional  space.  Parseval  s  theorem  is 

f  fh  <  i  v  )/  (.r  v  )dx  .dy  —  — f  Jh  P,  P,  )d  (3,  d P,.  (3  15) 

4tt" 

If  h<  vv  I  is  a  tomb  function  distributed  along  the  spatial  lattice  with  a  cell  to  cell  phase 


h  1  »  v  )  =  Z  Z  ^  P~P  .... 


(3  16a) 


i  hen  H  P  3  1  it  also  a  vornb  f  unction  distributed  along  the  reciprocal  lattice 


P  P  )  =  ~W-  Z  E  A<A-A'„r-< 


(  J  16b) 


3pplic.it  urn  >1  Parse-,  a  I  s  t  heorem  to  the  second  sum  of  1  uualion  t  1  1  4  )  yields 
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„  —jt  0  nA?o“?  1  2  +  (  I  z  u~:  ' 1  +cC-‘'  )2 

>)  =  — 4 - 

4ir  V  I  p{)-jy  1 2  +  (  1 2(I-J  '  I  +cCA  )2 

“7*  0  n/  1  PO_P  I  2  +  (  I  ;  0-z  '  I  )} 

- 1 - e~y  (f3'  *  <+'*  dxdy 

4ir  V  I  P,t-P  I  2  +  ( I  r  n-z  '  I +cCA  )2 

Substitution  of  Equations  (3.16)  and  (3.17)  into  (3.15)  yields 


m  a  >  =  // 


(3.17a) 


(3.17b) 


477“ 


Ioo  oo  -/ vlz-’l  4-cC  .7  — 

_  £  £  g  y  )(PP-P) 


r  A 

v *'*  '  ffl  =— oon  =— oo 


2y  v 


where  y  is  defined  by  Equation  (3.13). 


(3.18) 


3.5  Numerical  Implementation  of  the  Spatial  Domain  Acceleration 

As  with  the  strip  problem,  the  time  needed  to  calculate  the  convolution  of  a  two- 
dimensional  current  pulse  with  the  periodic  Green's  function  tested  with  a  delta  function 
was  plotted  for  three  different  methods.  In  this  section,  c  is  multiplied  by  the  unit  cell  area 
(C  A  ) 

Method  1  implements  the  integration  of  the  periodic  Green's  function  (Equation 
(3  14))  over  a  two-dimensional  patch  numerically  using  a  Romberg  integration  routine.  For 
every  x'.y'.z'  chosen  by  the  routine,  the  spectral  and  spatial  domains  of  Equation  (3.14) 
are  summed  to  accuracy  When  the  test  is  coincident  with  the  basis  function,  the  singular¬ 
ity  is  removed  from  the  non-asvmptoiic  term  and  added  back  in  analytically.  The  singu¬ 
larity  does  not  occur  in  the  asymptotic  terms  since  c>0  For  each  x  '.v  ‘.z  '.  the  test  point  of 
the  asymptotic  function  moves  as  shown  in  Figure  3  5  This  is  allowed  since  the  summa¬ 
tions  in  the  spectral  and  spatial  domains  are  being  done  together.  Figure  3  6  shows  the  time 
behavior  ol  Method  1  for  a  Imx  lm  basis  arranged  on  a  regular  hexagonal  lattice.  The 
time  required  to  calculate  the  matrix  element  for  a  plate  array  is  similar  to  the  time 
required  lor  a  strip  array  If  c  is  too  small  ( c  <0  02 ).  the  spearal  domain  is  overweighted, 
and  il  i.  is  too  large  (c  >0  12).  the  spatial  domain  is  overweighted  The  time  needed  to  calcu- 
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Figure  3.6  Time  required  for  Method  1  vs  c 
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late  the  matrix  element  for  the  plate  array  is  approximately  the  time  needed  to  calculate 
the  strip  array  matrix  element  raised  to  the  fourth  power. 

Method  2  involves  distributing  the  basis  and  testing  functions  onto  the  spatial  and 
spectral  domains. 


CO  QO 
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(3.19) 
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The  spatial  domain  integration  is  performed  numerically.  The  spectral  domain  integration 
is  performed  analytically.  When  there  is  no  overlap  in  z  between  J  and  T  .  the  spectral 
domain  contribution  becomes 
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With  overlap  the  spectral  domain  contribution  becomes 
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where 
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(3.20b) 


A  1  0,  COS0+/3*  sinO 

Since  the  spectral  and  spatial  domain  contributions  are  integrated  separately,  the  test  point 
of  the  asymptotic  terms  must  be  fixed  for  the  entire  calculation.  The  asymptotic  lest  point 
must  not  fall  on  the  basis  function  since  the  singularity  of  the  asymptotic  term  is  not  taken 
into  account.  Method  3  is  the  same  as  Method  2  with  a  point  approximation  for  the  spatial 
domain  integral  outside  the  core  region  (m-1 :1  ,n-l  .1 ). 

Figures  3.7-3.10  show  the  calculation  time  needed  for  all  three  methods  for  various 
locations  of  testing  functions.  In  all  cases.  Method  3  is  fastest  at  the  cost  of  accuracy 
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f  is;urr  3  9  Tim*  required  for  Methods  1.2  and  3  vs.  c  for  position  3  in  Figure  3.11 


Figure  3.1 1  Positions  of  testing  point  for  Figures  3  7-3 
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Me»h..d  1  is  most  accurate  lor  all  values  of  i  In  figure  3  7  the  situation  corresponds  to  1 
m  I  core  I  l  I  Method  I  shows  optimum  time  behavior  in  the  0.01  <t<0  1  range  where 

*  **  t  *;  h  t  mg  -!  speitra.  anil  spatial  domain  are  approximately  equal  Methods  2  and  3  show  a 

•  tr  p  n  me  jt  i -  ' ih  *  here  the  spectral  domain  test  point  moves  of!  plane  and  gets 

e\p.-nero  :jl  mn  .  ergense  figure  sorresponds  to  2  in  figure  3  II  The  test  point  is  "off 
plane'  t - 1  'yg  n  *itn  \s  a  result  son.ergence  tor  all  methods  is  more  rapid  than  in  Figure 

'  '  f  igure  \  -s;  he  sell  'erm  correspi>niling  to  1  in  I  igure  3  1  1  The  large  limes  here  arise 

1  ■  m  -riguiant  ies  ,p  t  he  as  .  mptot  n  terms  sime  until  i  exseeds  0  08.  the  testing  point  of  the 
i  -n  pi  '  ,i  1  e'rr>  t  a  I  Is  on  t  he  basts  t  uni  t  ion  I  his  san  be  seen  more  dramatically  if  Method 
'  •  - 1 a  rn  ned  a .  ne  I  gure  t  10) 
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4.  CONCLUSIONS 


This  report  has  investigated  the  convergence  characteristics  of  the  periodic  Green  s 
function.  The  characteristics  of  the  periodic  Green's  function  have  been  discussed  in  the 
framework  of  two  examples:  (i)  scattering  from  a  one-dimensional  array  of  strips,  and  (ii) 
scattering  from  a  two-dimensional  array  of  plates. 

The  periodic  problem  may  be  formulated  in  terms  of  responses  to  line/point  sources 
(spatial  domain)  or  in  terms  of  response  to  current  sheets  (spectral  domain).  The  spatial 
domain  is  slowly  convergent  everywhere  while  the  spectral  domain  is  only  slowly  conver¬ 
gent  in  the  "on  plane*  case  when  the  testing  function  is  located  in  the  array  plane.  The  slow 
convergence  in  one  domain  stems  from  an  unavoidable  singularity  in  the  reciprocal  domain. 

If  the  basis  function  is  located  entirely  in  the  array  plane  (flat  case),  the  derivatives 
may  be  transferred  onto  the  Green's  function  and  the  smoothness  of  the  basis  and  testing 
functions  may  be  used  to  help  convergence.  If.  on  the  other  hand,  the  basis  function  is 
rotated  in  the  array  plane,  then  all  derivatives  in  the  expression  must  be  transferred  expli¬ 
citly  onto  the  basis  and  testing  functions.  The  transfer  of  derivatives  in  the  scalar  potential 
term  and  approximations  in  the  vector  potential  term  simplify  the  problem  to  one  of 
finding  the  response  at  a  point  due  to  an  array  of  one-dimensional  current  pulses  in  the  case 
of  strips  and  two-dimensional  current  pulses  in  the  case  of  plates. 

In  order  to  quickly  do  the  summation  to  be  computed  in  TR  *  J  *  G  .  a  combination 
of  both  the  spectral  and  spatial  domains  must  be  used.  Accelerating  the  spatial  domain  is 
shown  symbolically  as 

T*  *  J  *  (G  — G°  )  +  F-,[G°  ]  (4.1) 

The  function  G  is  slowly  convergent  but  peaked.  A  smooth  function  that  asymptotically 

approaches  G  (G‘* )  is  subtracted  from  G  to  render  the  first  two  terms  in  the  brackets 
rapidly  converging.  The  function  G "  is  then  added  in  the  spectral  domain  using  the  Poisson 
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summation  formula.  Since  the  Poisson  summation  formula  essentially  Fourier  transforms 
a  smooth  function  with  wide  support,  this  term  is  also  rapidly  convergent  The  same  pro¬ 
cedure  may  be  applied  to  the  spectral  domain. 


■j~R  *  J  * 


F~KG-Ga  )  +  G" 


(4.2) 


The  function  C  in  Equation  (4.1)  was  chosen  by  moving  the  testing  point  "off  plane* 
through  use  of  a  parameter  c.  In  this  report,  c  is  multiplied  by  the  area  of  the  unit  cell  to 
get  a  distance.  Numerical  experiments  were  performed  to  determine  the  value  of  c  required 
to  minimize  the  time  needed  to  calculate  Equation  (4.1).  Three  methods  were  studied: 
Equation  (4.1)  itself  (Method  1):  distributing  the  basis  and  test  convolutions  onto  each 
domain  and  performing  the  convolution  analytically  in  the  spectral  domain  (Method  2); 
finally,  calculating  the  out  of  core  terms  of  the  spatial  portion  of  Method  2  using  a  point 
approximation  to  the  integrals  (Method  3). 


Method  1  is  the  most  accurate  of  all  the  methods  for  all  values  of  c  chosen.  The 
optimum  value  of  c  for  Method  1  is  in  the  range  0.01  < c  < 0  1.  For  this  range,  the  spatial 
and  spectral  domains  are  weighted  approximately  evenly.  Method  3  is  the  least  accurate  of 
the  methods  and  its  accuracy  decreases  as  c  is  increased,  due  to  the  approximation  in  the 
spatial  domain.  Method  3  is  also  the  fastest  method  of  the  three  for  a  wide  range  of  c.  The 
optimum  value  of  c  for  Method  3  due  to  the  accuracy  is  0.001  <c <0.03.  Method  2  has 
accuracy  problems  whenever  the  "on  plane*  case  occurs.  It  is  also  the  slowest  of  all  the 
methods  due  to  problems  in  specifying  the  absolute  accuracy  of  the  summations. 


In  summary.  Method  1  is  recommended  when  accuracy  is  the  prime  concern  while 
Method  3  is  recommended  when  speed  is  desired.  In  all  cases,  the  choice  of  c  must  be  made 
to  ensure  that  the  asymptotic  term  test  point  does  not  fall  on  the  basis  function  since  the 
singularity  of  the  asymptotic  term  has  not  been  taken  into  account. 


In  general,  it  was  found  that  when  the  strips  or  plates  are  flat  in  the  array  plane,  the 
smoothness  and  width  of  the  basis  and  testing  functions  help  the  convergence  of  the 


spectral  domain  so  much  that  the  spectral  domain  should  get  the  entire  weighting 


Acceleration  techniques  need  not  be  applied  In  the  cases  where  the  strips  or  plates  are 
rotated  with  respect  to  the  array  plane  acceleration  techniques  can  be  applied  which  results 
in  a  substantial  time  saving  when  c  is  selected  in  the  ranges  recommended  above 
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